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We study evolutionary game dynamics on structured populations in which individuals take part in several 
layers of networks of interactions simultaneously. This multiplex of interdependent networks accounts for 
the different kind of social ties each individual has. By coupling the evolutionary dynamics of a Prisoner's 
Dilemma game in each of the networks, we show that the resilience of cooperative behaviors for extremely 
large values of the temptation to defect is enhanced by the multiplex structure. Furthermore, this resilience 
is intrinsically related to a non-trivial organization of cooperation across the network layers, thus providing 
a new way out for cooperation to survive in structured populations. 

The understanding of the emergence of cooperative behavior in human and animal societies 1 ' 2 as well as in 
other contexts (e.g., the formation of multicellular organisms or their organs 3 ) is a major challenge in 
science. Interdisciplinary physicists have paid attention to this problem because of the underlying nonlinear 
and stochastic nature of the interactions among the entities involved 4-6 . The mathematical setting that has led to 
many deep insights about this problem is evolutionary game theory 7 " 9 , that allows to formulate in quantitative 
terms the most important prototypical social interactions, such as conflicts and/or dilemmas 10 . Scientists have 
unraveled that a key issue to ascertain the evolution of cooperation is the network of relationships 11 " 13 among the 
intervening agents. This drive us to evolutionary game theory on graphs, one of the most intriguing dynamical 
processes on networks and one that is currently receiving a lot of attention 14 " 21 . 

As network science evolves 22 " 25 , new questions about the capital problem of the emergence of cooperation arise. 
The empirical resolution of the structure and time evolution of social ties has been simultaneously improved. 
These advances has been largely facilitated by the explosion of data about mobile communication 28 29 , web -based 
social platforms 26 ' 27 and even the monitoring of face-to-face human interactions 30 ' 31 . Thus, although the network 
perspective has offered a novel way out for cooperation to survive in social systems 32 " 35 , the latter advances on the 
characterization of social systems demand more work to unveil the influence that social patterns have on the 
evolution of cooperation. Particularly important in the description of social systems are those structures that 
account for multiple types of links and time-evolution of links, commonly known as multiplex 36 " 38 . 

Social systems are shown as a superposition or projection of a number of interdependent complex social 
networks, where nodes represent individuals and links account for different kind of social ties such as those 
stablished with family relatives, friends, work collaborators, etc. In our daily life we experience this social splitting 
by distinguishing our behavior within each of the social layers we belong to. However, the influence that the 
multiplex nature of social interactions has on the evolution of cooperation is still an open question, being recently 
tackled 39 within a framework consisting in two coupled networks. Besides, on more general grounds, it has been 
recently shown that the interdependent structure can influence dramatically the functioning of complex systems 
in the context of percolation 40 " 44 and cascade failures 41 ' 45 . Thus, it is necessary to study how the interplay among 
such multiple interdependent social networks affects the onset of large scale human behavior, in particular the 
emergence of cooperation in multiplexes. 

In this report, we address the problem of the emergence of cooperation in multiplexes. We will use tools of 
complex networks and evolutionary game theory to shed light on the emergence of cooperation in populations of 
individuals participating simultaneously in several networks in which a Prisoner's Dilemma game is played. Our 
results show that a mutiplex structure enhances the resilience of cooperation to defection. Moreover, we show that 
this latter enhancement relies on a nontrivial organization of the cooperative behavior across the network layers. 

Let us first describe the multiplex backbone in which the evolutionary game is implemented. We consider a set 
of m interdependent networks each of them containing the same number AT nodes and L links. Each individual is 
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Figure 1 | Schematic representation of a multiplex network. The 

multiplex is made of N = 5 nodes embedded within m = 3 interdependent 
networks (or layers) each one containing L = 3 links. 

represented by one node in each of the m networks layers while its 
neighbors are (in principle) different for each of the layers (see Fig. 1). 
In this way we define a set of m adjacency matrices {A 1 } (with 1=1, 
. . ., m) so that A\j = 1 when individuals i and; are connected within 
network / whereas A ■■ = 0 otherwise. Thus, a given individual, say i, is 
connected to k\ = ^ A\- other individuals within network layer /. In 
our case we will consider that each of the layers is an Erdos-Renyi 
(ER) random graph characterized by an average degree (k) = 2L/N. 
In this way, the probability that an individual is connected to k 
individual in a given layer is given by the Poisson distribution: 
P(k) = (k) k exp(-fc)/fc. 

Having introduced the multiplex composed by the set of m inter- 
dependent networks we now focus on the formulation of the evolu- 
tionary dynamics. Each of the individuals, say i, adopts a given 
strategy s\{t) for playing with its neighbors in network / at time step 
t. This strategy can be cooperation \s\{t) = l] or defection [s\(t) = 0] . 
Then, at each time step, each individual plays a Prisoner's Dilemma 
(PD) game with its neighbors in network /. For each of the k\ PD 
games played within network layer / an individual i facing a co- 
operator neighbor will collect a payoff 1 or b > 1 when playing as 
cooperator or as defector respectively. On the contrary if i faces a 
defector it will not collect any payoff regardless of its strategy. This is 
the weak version of the PD game which makes use of a single para- 
meter b accounting for the temptation of playing as defector. After 
round t of the PD game, an individual has played once with its k\ 
neighbors in layer / thus collecting an overall payoffpj(f). Obviously, 
the net payoff of a player i is the sum of all the payoffs collected in 
each of the m network layers, Pj(t) = p\> achieved by using a set of 
strategies {sj(f)}- 

Once the PD is played, all the players update their strategies 
simultaneously, i.e., we consider synchronous updates. The update 
process makes use of the replicator-like rule that works as follows. 
Each of the players, say i, chooses a layer, say /, at random among the 
m possible networks and a neighbor j (also randomly) among its k\ 



acquaintances. Then it compares their total payoffs, P,-(f) and Pj(t), 
obtained in the last round of the game. If P,-(f) > P ; (£) nothing 
happens and i will use the same strategy within the network layer / 
in the next round of the PD game, Sf(t + 1) = s,-(f). However, when 
Pj{t) > Pi(t) agent i will take the strategy of j at layer / with a 
probability proportional to their payoff difference: 



b max(Ki,Kj) ' 



(1) 



where K { = J2i^\- Note that the update process entangles the evolu- 
tionary dynamics of the network layers as the choice of what strategy 
will be used in a given layer during the next round of the game 
depends on the overall payoffs, not only on the payoffs obtained 
in the particular layer. Although other entanglements are also pos- 
sible to make the evolutionary dynamics of the m networks inter- 
dependent, the one used in this work relies on the social nature of 
layers' interdependency. While two neighbors within a given layer 
know the strategies used by each other in the layer, they are unaware 
of their opponent's strategies in the remaining networks, so they 
have to assume that the total benefit or success achieved by their 
layer's acquaintances is the outcome of using the observed strategy in 
each of the m networks of the multiplex. 

Following the above evolutionary rules we let evolve the states of each 
individual in each of the layers and compute the instantaneous level of 
cooperation d{t) in each layer / and in the whole multiplex c(t) as: 



, m -, m N 



(2) 



Results 

In Fig. 2 we plot the average fraction of cooperators (c) (see Methods) 
versus b for different values of the number m of layers and two values 
of the layers' average degree (k) = 3 and 20. Note that the case m = 1 
corresponds to the absence of layers' interdependency. While for low 
values of the temptation b the average level of cooperation on the 
multiplex decreases with the number m of layers, it increases with m 
for higher values of b, so that the decrease of cooperation with b 
becomes progressively slower as the number of layers in the multi- 
plex increases. Importantly, the resilience of cooperation observed 
for large values of the temptation to defect is not restricted to the 
weak version of the PD game; we have checked that our results 
qualitatively remain unaltered when two defectors playing get a pos- 
itive payoff, i.e., when defection is a strict best response to itself. 

The possibility that agents adopt different strategies in different 
layers is crucial for the resilience of cooperation in multiplexes, as 
revealed in Fig. 3, where we show the average cooperation level 
achieved in the multiplex when agents use homogeneous strategies 
(s- independent of /), i.e., the same strategy in the m layers. Indeed, for 
homogeneous strategists' populations the cooperation decays dra- 
matically with b, the faster the larger the value of m. This can be easily 
understood as follows: The situation when agents use homogeneous 
strategies is equivalent to consider the standard evolutionary game 
dynamics on a network whose set of links is the union of the set of 
layer's links (with weights assigned to them whenever two agents are 
neighbors in more than one layer in the multiplex). This network has 
(approximately) an average degree of m(k), which for large values of 
m approaches a well-mixed population, where cooperation extin- 
guishes quickly. 

The survival of cooperation in multiplexes is supported by the 
network structure inside the multiplex layers, i.e. on network reci- 
procity 46 . One can easily prove that if the layers are assumed to be 
well-mixed (fully connected), then the only surviving strategy of the 
evolutionary dynamics is to defect in all the layers. To characterize 
the degree of heterogeneity of the surviving strategies on the 
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Figure 2 | Cooperation diagrams of multinetworks. Average level of cooperation (c) as a function of the temptation to defect b for several multinetworks 
with different number of layers m. In panel A the network layers are ER graphs with (k) = 3 (sparse graphs) while in panel B we have (k) = 20. In both cases 
N = 250 nodes. As can be observed, the resilience of cooperation increases remarkably as the number of layers m grows. Finally, panel C shows the curves 
(c)(b) for ER graphs with (k) = 3 (as in panel A) for m = 2 and m = 10 and different network sizes N = 100, 200 and 400. 



multiplex, let us consider the fraction Xj(t) of layers where agent i 
plays as a cooperator at time t: 



I m 



(3) 



that, after averaging over the observation time interval T (see above) 
defines the variable x, whose probability density P(x) as a function of 
b is shown in Fig. 4. One should remark that P(x) exhibits a well 
defined maximum for each value of b, along with a relatively small 
width around it. In other words, for each value of b there is a truly 
characteristic value of the fraction x of layers where a randomly 
chosen agent behaves cooperatively. On the contrary, the layer's 
cooperation level c l possess a rather wide distribution density, as 
shown in Fig 5. In a given multiplex realization, and for a given 
value of the temptation parameter b, layers with quite dispersed 



cooperation levels coexist, so that there is no truly characteristic value 
for the layer's cooperation. In other words, the cooperative compo- 
nents [whose number is typically m-x max , where x max is the location of 
the maximum of P(x)] of the agents' set of strategies {s\(t) } are not 
uniformly distributed on the multiplex layers. 

The above findings suggest the failure of "mean-field-like" theor- 
etical explanations of the observed behavior based on the absence of 
correlations. As an example, consider a layer /, and let us assume that 
the payoff an agent receives from the games played on the rest of 
layers is normally distributed according to the hypothesis that it has 
(k) neighbors on each layer and that agents use the cooperative 
strategy in each layer with a probability x, independently of each 
other. One can compute easily the mean and the variance of this 
payoff distribution, as a function of x, m and b (see Methods). For 
self- consistency, one can assign to x the value c\t) of the instantaneous 
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Figure 3 | Comparison with the cooperation diagrams of null models. Average degree of cooperation (c) as a function of b (solid line with filled circles) 
for m = 3 (A), 6 (B), 10 (C) and 20 (D). In each panel we show the case of a simplex (m = 1) network, the evolution (c)(b) for homogeneous strategists' 
populations corresponding to each value of m (Homo.) and the curve (c)(b) corresponding to the mean-field assumption for the coupling between layers 
(M.F.). Each multiplex network has N = 250 nodes while the interdependent layers are ER graphs with (k) = 3. 
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cooperation in the considered layer, and then proceed to run the 
evolutionary game dynamics on the layer with payoffs given by the 
sum of the payoff obtained in the layer plus a payoff taken from this 
"mean -field-like" distribution. The results obtained from this 
approach are shown in Fig. 3. The only feature of the behavior 
observed in the multiplex that seems to be slightly captured, but only 
qualitatively, by this approach is the tendency when m increases to the 
decrease of cooperation for low values of b, and its increase for high 
values of b. The building up of correlations in the distribution over the 
multiplex layers of the fraction of cooperative strategies appears to be 
an essential ingredient for the observed enhancement of cooperation 
in multiplexes. 
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Figure 4 | Evolution of the degree of cooperation of individuals across 
layers. Each contour plot shows the evolution of the probability P(x) of 
finding an individual playing as cooperator in a fraction x of the network 
layers as b increases. In A m = 15 while for B m = 20. In both cases the 
networks have N = 250 nodes while each layer is an ER graph with (k) = 3. 
Observe that the resilience of cooperation is intrinsically due to the fact 
that individuals play different strategies across the different layers (no 
homogeneous strategists appear until defection dominates at very high 
values of b). For each value of b there is a well defined maximum for P(x). 
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Discussion 

Summing up, we have incorporated the multiplex character of social 
interactions into the formulation of evolutionary games in structured 
populations. By considering a Prisoner's Dilemma game we have 
shown that cooperation is able to resist under extremely adverse 
conditions, for which the usual simplex formulation, i.e. the network 
approach, fails. In a nutshell, the addition of network layers has two 
effects. On one hand the level of cooperation for low values of the 
temptation to defect appears to decrease. However, the enhancement 
of the cooperation resilience shows up when temptation is further 
increased. It is in the region of large temptation when the inter- 
dependency of network layers outperform, regarding the average 
level of cooperation, the behavior found in simplex networks reach- 
ing values of b > 3 for which even scale-free networks fail to sustain 
cooperation 34,35 . 

The observed resilience of multiplex networks is sustained in the 
segregation of cooperative and defective strategies across the mul- 
tiple network layers contained in the multiplex. Moreover, we have 
shown that this segregation is non -trivial by comparing with mean- 
field approaches producing no cooperation. Thus, our results point 
out a complementary mechanism to the so-called network recipro- 
city, paving the way to the study of more complex and realistic 
multiplex architectures and alternative dynamical couplings between 
the networks embedded in them. Let us note, that we have made use 
of network layers with a regular (Erdos-Renyi) topology in order to 
avoid spurious effects, such as the degree heterogeneity, that may 
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Figure 5 | Histograms of layers' cooperation level. The panels show the probability of finding a level of cooperation c l in a randomly chosen layer of the 
multiplex. The results are obtained from a multiplex of m = 20 layers and four different values of the temptation parameter b= 1.0, 1.2, 1.4, and 1.6 from 
panels A to D. 50 realizations of the multiplex were employed. 
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contribute to the enhancement of cooperation. However, the study of 
multiplex topologies incorporating the interdependency of scale-free 
layers seems a promising continuation of the results presented here. 

Methods 

Numerical and statistical details of the simulations. In our simulations we start 
from a configuration in which a player i in layer / cooperates or defects with equal 
probability. Then, we run the evolutionary dynamics for a transient time t 0 of typically 
t 0 = 2X 10 4 generations. After this transient period we further iterate the evolutionary 
dynamics over a time interval T of typically 10 5 generations. It is during this latter 
window when we compute the quantities of interest such as the average cooperation 
level in a given layer /, c l (l = 1, . . . ,m): 

t=t 0 i=l 

or the average cooperation in the whole multiplex c: 

1 m 

The values shown in each plot represent the average of the above quantities over a 
number of realizations, typically 50. After this latter average we obtain the final 
average level of cooperation (c) for each value of the temptation parameter b. 

Mean-field calculation. The mean-field assumption for the coupling of layers 
assumes that there is no correlation between the strategies used by an individual in 
each of the m layers. To this aim, we consider a network (single layer) in which an 
individual plays a PD game with his neighbors receiving a payoff. In addition we add 
to this latter payoff a quantity mimicking the payoff obtained in the rest of (m — 1) 
layers. This additional payoff is randomly assigned after each round of the game from 
a normal distribution whose precise form depends on the number of cooperators in 
the system. 

To compute the mean and variance of the normal distribution at work, we first 
consider that all the nodes in the network are connected to (k) neighbors in each of the 
network layers. In this way, the possible payoffs of a given individual i in one of the 
network layers /are p\ = {0,1,2, . . . ,(k),b,2b, . . . ,(k)b}. Considering now that players 
use the cooperative strategy in each layer with probability x, independently of each 
other, and that the multiplex is composed of m network layers we can assign compute 
the probability q (p{) that player i obtain a payoff p\ in layer /. In our case, we fix (k) = 3 
so that there are 7 possible payoffs whose probabilities read: 
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= x(}-xf + (}-x)\ 
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= 3(xf(l-xf, 
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-3b) 


= (xf(\-x). 



With these expressions one can easily compute the expected value for the payoff 
obtained by an individual i in a layer / given the value of x as: 

7 

phE^)^- 3 ^ 1 ( 6 ) 

while the variance of the above expected value reads: 

^=E^(^)^-^ = 3^{(l + 2x)[^ 2 (l-x)+x]-3xWl-x)+x] 2 }. (7) 
j=i 

Given the values of x and b, equations (6) and (7) allow us to compute the normal 
distribution and assign the additional payoffs a player receives from the other 
(m — 1) layers. For self- consistency, at each time step, the value of x is re-computed 
from the fraction of cooperators in the system, i.e. c\t), in the considered layer. 
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